﻿# encodings UTF-8
import numpy.fft as nf
import numpy as np
import matplotlib.pyplot as plt
import os
import linecache
import re

plt.rcParams["font.sans-serif"] = ['SimHei']
plt.rcParams["axes.unicode_minus"] = False

path = 'Data'  # 数据存储的文件夹名
files = os.listdir(path)
for file in files:
    if file.startswith(".") or file.endswith("1") or file.endswith(".png") \
            or file.endswith(".pdf") or file.endswith(".txt"):  # 过滤掉以“.”开头的文件和以“1”、“.png”结尾的文件
        continue

    # scale = float(re.split(r'[\t\n]', linecache.getline(path + r'/' + file, 2))[1])  # 读取电压信号的比例值
    scale = 1
    data = np.loadtxt(path + r'/' + file, skiprows=9)  # 跳过前9行数据
    t = data[:, 0]  # 第一列数据为时间序列
    y = data[:, 1]  # 第二列数据为函数序列

    '''
    对自己定义的函数进行分析
    '''
    # B = 30.0    # 频谱分析的最大频率
    # f_s = 2 * B # 采样频率
    # delta_f = 0.01  # 频谱分辨率
    # N = int(f_s / delta_f)  # 采样点数
    # T = N / f_s # 采样时间

    # def signal_samples(t):
    #     # return (2 * np.sin(2*np.pi*t) + 3 * np.sin(22*2*np.pi*t) + 2 * np.random.randn(*np.shape(t)))
    #     return 2 * np.sin(4*2*np.pi*t) + 3 * np.sin(3*2*np.pi*t)

    # t = np.linspace(0, T, N)
    # f_t = signal_samples(t)

    fig, axes = plt.subplots(1, 2)  # 画两个子图，分别为时域图和频域图
    axes[0].plot(t, y)
    axes[0].set_xlabel("time (s)")
    axes[0].set_ylabel("signal")

    Y = nf.fft(y)
    N = y.size  # 采样点数
    delta_t = t[1] - t[0]  # 采样间隔
    f = nf.fftfreq(y.size, delta_t)  # 频域的横坐标序列
    mask = np.where(f > 0)  # 只分析频率大于零的部分（小于零的部分与之完全对称）
    A = scale * 2 * abs(Y[mask]) / N  # 频谱图的纵坐标变换为实际振幅

    meas_freq = float(file.split('Hz')[0])  # 从文件名中读取激振器的频率值
    sign_freq = (f[f > 0])[np.argmax(abs(Y[mask]))]  # 找到频谱图山振幅最大时对应的频率值
    amplitude = np.max(A)  # 最大振幅值
    print("\n分量最大的频率分量为：{:.2f} Hz, 对应的幅值为：{:.2f} V".format(sign_freq, amplitude))
    devi_freq = sign_freq - meas_freq  # 绝对误差
    real_devi = devi_freq / meas_freq  # 相对误差
    print("偏差为：{:.2f} Hz（{:.2%}）".format(devi_freq, real_devi))
    text = file + "\n" + 'f = {:.2f} Hz, A = {:.2f} V, error = {:.2f} Hz, {:.2%}'.format(sign_freq, amplitude,
                                                                                         devi_freq, real_devi)
    plt.suptitle(text)  # 将误差作为图的标题，显示在图的上方

    axes[1].plot(f[mask], A)
    axes[1].set_xlabel("frequency (Hz)")
    axes[1].set_ylabel("$A$")
    plt.tight_layout(rect=[0, 0, 1, 1])  # 设置子图布局，避免重叠
    plt.savefig(path + r'/' + file + '.pdf')  # 将分析结果保存在以频率信息命名的png图片中

plt.show()
